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Abstract. A finite volume scheme for the (Patlak-) Keller-Segel model in two space di- 
mensions with an additional cross-diffusion term in the elliptic equation for the chemical 
signal is analyzed. The main feature of the model is that there exists a new entropy func- 

£S) tional yielding gradient estimates for the cell density and chemical concentration. The 

,_i main features of the numerical scheme are positivity preservation, mass conservation, en- 

tropy stability, and — under additional assumptions — entropy dissipation. The existence 
of a discrete solution and its numerical convergence to the continuous solution is proved. 
Furthermore, temporal decay rates for convergence of the discrete solution to the homo- 
geneous steady state is shown using a new discrete logarithmic Sobolev inequality. Nu- 

i— i merical examples point out that the solutions exhibit intermediate states and that there 

exist nonhomogeneous stationary solutions with a finite cell density peak at the domain 

^ | boundary. 

-i— > 

a 

1. Introduction 

Chemotaxis, the directed movement of cells in response to chemical gradients, plays 
an important role in many biological fields, such as embryogenesis, immunology, cancer 
growth, and wound healing [2T1 I3T] . At the macroscopic level, chemotaxis models can 
be formulated in terms of the cell density n(x,t) and the concentration of the chemical 
signal S(x,t). A classical model to describe the time evolution of these two variables is 
the (Patlak-) Keller-Segel system, suggested by Patlak in 1953 [29] and Keller and Segel 
in 1970 [23]. Assuming that the time scale of the chemical signal is much larger than that 
of the cell movement, the classical parabolic-elliptic Keller-Segel equations read as follows: 

> d t n = div(Vn - nVS), = AS + fin - S in Q, 

where Q C IR 2 is a bounded domain or Q = IR 2 , with homogeneous Neumann boundary 
and initial conditions. The parameter /i > is the secretion rate at which the chemical 
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substance is emitted by the cells. The nonlinear term riVS* models the cell movement 
towards higher concentrations of the chemical signal. 

This model exhibits the phenomenon of cell aggregation. The more cells are aggregated, 
the more the attracting chemical signal is produced by the cells. This process is counter- 
balanced by cell diffusion, but if the cell density is sufficiently large, the nonlocal chemical 
interaction dominates diffusion and results in a blow-up of the cell density. In two space 
dimensions, the critical threshold for blow-up is given by M = j n n dx = 4ir if Q is a 
bounded connected domain with C 2 boundary [27] and M = 87r in the radial and whole- 
space case [3],|28]. The existence and uniqueness of smooth global-in-time solutions in the 
subcritical case is proved for bounded domains in [23] and in the whole space in [SJ. In 
the critical case M = 8tt, a global whole-space solution exists, which becomes unbounded 
as t — > oo [1]. Furthermore, there exist radially symmetric initial data such that, in the 
supercritical case, the solution forms a (^-singularity in finite time [20] . 

Motivated by numerical and modeling issues, the question how blow up can be avoided 
has been investigated intensively the last years. It has been suggested to modify the 
chemotactic sensitivity (modeling, e.g., volume-filling effects), to allow for degenerate cell 
diffusion, or to include suitable growth-death terms. We refer to [22] for references. Another 
idea is to introduce additional cell diffusion in the equation for the chemical concentration 
[HI 122] • This diffusion term avoids, even for arbitrarily small diffusion constants, the blow- 
up and leads to the global-in-time existence of weak solutions [22]- The model, which is 
investigated in this paper, reads as follows: 



where 5 > is the additional diffusion constant. We impose the homogeneous Neumann 
boundary and initial conditions 



The advantage of the additional diffusion term is that blow-up of solutions is translated 
to large gradients which may help to determine the blow-up time numerically. Another 
advantage is that the enlarged system ([I]) exhibits an interesting entropy structure (see 
below) . 

At first sight, the additional term 6 An seems to complicate the mathematical analysis. 
Indeed, the resulting diffusion matrix of the system is neither symmetric nor positive 
definite, and we cannot apply the maximum principle to the equation for the chemical 
signal anymore. It was shown in |22j that all these difficulties can be resolved by the 
observation that the above system possesses a logarithmic entropy, 



(1) 



d t n = div(Vn - nVS), = AS + 5 An + fin- S, x G Q, t > 0, 



(2) 



Vn ■ v = VS • v — on dQ, t > 0, n(-, 0) = n in Q. 



E(t) 




which is dissipated according to 




A\Vjn\ 2 + \\^S\ 2 + \s 2 \dx 



5' 1 8 J 



s 



L 

Jn 



nSdx. 
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Suitable Gagliardo-Nirenberg inequalities applied to the right-hand side lead to gradient 
estimates for y/n and S, which are the starting point for the global existence and long-time 
analysis. 

In this paper, we aim at developing a finite volume scheme which preserves the entropy 
structure on the discrete level by generalizing the scheme proposed in [16]. In contrast to 
[16] , we are able to prove the existence of discrete solutions and their numerical convergence 
to the continuous solution for all values of the initial mass. Moreover, we show that the 
discrete solution converges for large times to the homogeneous steady state if /i or 1/5 are 
sufficiently small, using a new discrete logarithmic Sobolev inequality (Proposition [3]). 

In the literature, there exist several approaches to solve the classical Keller- Segel sys- 
tem numerically. The parabolic-elliptic model was approximated by using finite difference 
[341 136] or finite element methods [25j [321 EH]. Also a dynamic moving-mesh method [6], 
a variational steepest descent approximation scheme [2] , and a stochastic particle approxi- 
mation [THl [19] were developed. Concerning numerical schemes for the parabolic-parabolic 
model (in which d t S is added to the second equation in 0), we mention the second-order 
central- upwind finite volume method of [11], the discontinuous finite element approach of 
[H], and the conservative finite element scheme of [33]. We also cite the paper [7] for a 
mixed finite element discretization of a Keller-Segel model with nonlinear diffusion. 

There are only a few works in which a numerical analysis of the scheme was performed. 
Filbet proved the existence and numerical convergence of finite volume solutions [16] . Error 
estimates for a conservative finite element approximation were shown by Saito [321 133] . 
Epshteyn and Izmirlioglu proved error estimates for a fully discrete discontinuous finite 
element method [14J. Convergence proofs for other schemes can be found in, e.g., [2| fl9]. 

This paper contains the first numerical analysis for the Keller-Segel model ((!]) with ad- 
ditional cross-diffusion. Its originality comes from the fact that we "translate" all the 
analytical properties of [22J on a discrete level, namely positivity preservation, mass con- 
servation, entropy stability, and entropy dissipation (under additional hypotheses). 

The paper is organized as follows. Section [2] is devoted to the description of the finite 
volume scheme and the formulation of the main results. The existence of a discrete solu- 
tion is shown in Section [3] A discrete version of the entropy-dissipation relation ^ and 
corresponding gradient estimates are proved in Section [4} These estimates allow us to 
obtain in Section [5] the convergence of the discrete solution to the continuous one when 
the approximation parameters tend to zero. A proof of the discrete logarithmic Sobolev 
inequality is given in Section [6] The long-time behavior of the discrete solution is investi- 
gated in Section [7j Finally, we present some numerical examples in Section [8] and compare 
the discrete solutions to our model ([TJ with those computed from the classical Keller-Segel 
system. 

2. Numerical scheme and main results 

In this section, we introduce the finite volume scheme and present our main results. 

2.1. Notations and assumptions. Let Q C 1R 2 be an open, bounded, polygonal subset. 
An admissible mesh of Q is given by a family T of control volumes (open and convex 
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polygons), a family S of edges, and a family of points (x^KeT which satisfy Definition 
9.1 in [13]. This definition implies that the straight line between two neighboring centers 
of cells (xk,%l) is orthogonal to the edge a = K\L. For instance, Voronoi meshes are 
admissible meshes [131 Example 9.2]. Triangular meshes satisfy the admissibility condition 
if all angles of the triangles are smaller than 7r/2 [T3J Example 9.1]. 

We distinguish the interior edges a G £ int and the boundary edges a G £ ex t- The set of 
edges £ equals the union £ int U £ ext . For a control volume K G T, we denote by 8k the 
set of its edges, by S intj K the set of its interior edges, and by £ e xt,A" the set of edges of K 
included in dfl. 

Furthermore, we denote by d the distance in IR 2 and by m the Lebesgue measure in IR 2 
or IR. We assume that the family of meshes satisfies the following regularity requirement: 
there exists £ > such that for all K G T and all a G &\t&,k with a = K\L, it holds 

(4) d(x K , cr)>£ d(x K , x L ). 

This hypothesis is needed to apply discrete Sobolev-type inequalities [1|. Introducing for 
a G E the notation 

d(x K , x L ) if a G £ int , a = K\L, 
d(x K ,a) if a G 8 ext>K , 



d a = 

we define the transmissibility coefficient 



m(cr) 

a G t. 



The size of the mesh is denoted by 

Ax = maxdiamffO. 

Let T > be some final time and Mt the number of time steps. Then the time step size 
and the time points are given by, respectively, 

T 

At = , t k = kAt, < k < Mt- 

We denote by T> an admissible space-time discretization of VLt = Q x (0,T) composed 
of an admissible mesh T of Q and the values At and Mt- The size of this space-time 
discretization T> is defined by rj = max{Ax, At}. 

Let X(T) be the linear space of functions Q — > IR which are constant on each cell K G T. 
We define on X(T) the discrete L p norm, discrete W 1,p seminorm, and discrete W 1,p norm 
by, respectively, 

i/p 




\ u \\o,p,t = ^ m(K)\u 



p 



V^tnt / 
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IM|l,P,T — IMkp.T + Ml,p,T ? 

where u G X(T), 1 < p < oo, and D a u = \uk — ul\ for a = K\L G £ in t- 

2.2. Finite volume scheme and main results. We are now in the position to define 
the finite volume discretization of (JT])-((2|) . Let D be a finite volume discretization of Qt- 
The initial datum n is approximated by its L 2 projection on control volumes: 

(5) n ^=y2 n °K 1 K, where n° K = 1 / n (x)dx, 

and 1^ is the characteristic function on K. Denoting by n k K and S K approximations of 
the mean value of n(-,t h ) and S(-,t k ) on K, respectively, the numerical scheme reads as 
follows: 

(6) m(JO^^* - E + E ^((DS^n^ - (DS™)- n ™) = 0, 

<76£a' ""Spirit, 

(7) - E - 5 E r Dn™ = m(K)(^ - ^ +1 ), 

for all .fT G T and < k < Mt — 1. Here, t> + = max{0, f }, t> ~ = max{0, — v}, and 

(8) DU k fUZ-U$ for a = K\LeS mt , K , 

' i for cr G t cxt ,x. 

The approximation S K is computed from ^ with /c = — 1. This scheme is based on a fully 
implicit Euler discretization in time and a finite volume approach for the volume variable. 
The implicit scheme allows us to establish discrete entropy-dissipation estimates which 
would not be possible with an explicit scheme. This approximation is similar to that in 
[16] except the additional cross-diffusion term in the second equation. 
The numerical approximations nx> and S*r> of n and S are defined by 



n v {x,t) = E n^ l l K (x), S v {x,t) = E S^ +1 l K {x), where x G Q, t G (t fe ,t A 



/.fc+ii 

A'er ifeT 

and = 0, . . . , Mt — 1. Furthermore, we define approximations V^n© and V^Sd of the 
gradients of n and S, respectively. To this end, we introduce a dual mesh: for K G T and 
a G 8k, let Ta: )0 - be defined by: 

• If a = K\L G £int,K, T Ka is the cell ("diamond") whose vertices are given by xk, 
Xl, and the end points of the edge a = K\L. 

• If a G £ext,K, T K>a is the cell ("triangle") whose vertices are given by xk and the 
end points of the edge a = K\L. 

An example of construction of Tx, a can be found in [TU]. Clearly, Tx a defines a partition 
of Q. The approximate gradient V^hd is a piecewise constant function, defined in = 
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Q x (0, T) by 

V v n v (x,t) = -^l-Dn k +^ K , a , x G T K , a , t G (t k ,t k+1 ), 

m{-LK,cr) 

where Dn k j^l is given as in ^ and UK,a is the unit vector normal to a and outward to K. 
The approximate gradient V v Sd is defined in a similar way. 

Our first result is the existence of solutions to the finite volume scheme. 

Theorem 1 (Existence of finite volume solutions). Let Q C M 2 be an open, bounded, poly- 
gonal subset and let T> be an admissible discretization of Q x (0,T). The initial datum 
satisfies n G L 2 (Vt), n > in Q. Then there exists a solution {(n^, S^-), K G T, < 
k < M r } to @-(|7|) satisfying 

(9) n k K > for all K G T, k> 0, 

(10) Yl < K >k = £ m (*K = 11^11^(0) forallk>0. 
KeT K&T 



Properties (|Qj) and ( 10 ) show that the scheme is positivity preserving and mass conser- 
ving. It is also entropy stable; see (17) below. 

Let (V v ) v>0 be a sequence of admissible space-time discretizations indexed by the size 
rj = max{ Ax, At} of the discretization. We denote by (7^)^>o the corresponding meshes 
of fi. We suppose that these discretizations satisfy Q uniformly in r], i.e., £ > does 
not depend on 77. Let (n v , S^) := (n^ , S"^ ) be a finite volume solution, constructed in 
Theorem [lj on the discretization D^. We set := V v,] . Our second result concerns the 
convergence of (n v , S v ) to a weak solution (n, S) to ([TJ-Q. 

Theorem 2 (Convergence of the finite volume solutions). Let the assumptions of Theorem 
[7] hold. Furthermore, let (D r? ) r?>0 be a sequence of admissible discretizations satisfying Q 
uniformly in rj, and let (n v , S^) be a sequence of finite volume solutions to (Js]) -([7]) . Then 
there exists (n, S) such that, up to a subsequence, 

71,7 — >■ n strongly in L 2 (f2<r), 

V^n^ — Vn weakly in L 2 (Qt), 

S v — S", ViS^ — VS weakly in L 2 (Q T ), 

and (n,S) G L 2 (0,T; H 1 ^)) 2 is a weak solution to in the sense of 

(11) / f{nd t (j)- Vn- V<f) + nVS ■V<f>)dx+ [ n (f>{-,0)dx = 0, 



(12) / I {-VS ■ V0 - SVn ■ V0 + nn<j> - S<p)dx = 

Jo Jq 

for all test functions (J) G C£°(fi x [0,T]). 

It is shown in J22J Theorem 1.3] that, if the secretion rate u > is sufficiently small or 
the diffusion parameter S > is sufficiently large, the solution (n, S) to @-(|2]) converges 



A FINITE VOLUME SCHEME FOR A KELLER-SEGEL MODEL 7 

exponentially fast to the homogeneous steady state (n*,S*), where n* = \\n \\ L i^ n )/m(Q) 
and S* = fin*. The proof in |22j is based on the logarithmic Sobolev inequality Therefore, 
we state first a novel discrete version of this inequality which is proved in Section [6j 

Proposition 3 (Discrete logarithmic Sobolev inequality). Let Q C M. d (d > 1) be an open 
bounded polyhedral domain and let T be an admissible mesh ofVL satisfying (|4]). Then there 
exists a constant Cl > only depending on Q, d, and £ such that for all u G X(T), 



f u 2 

/ ^ 2 log — _ 1N ||2 — dx < Cl\u\1 2T , 
Jn m IMIo,2.r 



where we abbreviated m = m(Q). 

The constant Cl > can be made more precise. Let Cs{q) > be the constant in the 
discrete Sobolev inequality [1, Theorem 4] 

(13) IM|o l9l r < ^^IM|i,2,t for u E X(T), 

where 1 < q < 2d/(d — 2) (and 1 < q < oo if d < 2), and let Cp(r) > be the constant in 
the discrete Poincare-Wirtinger inequality [1, Prop. 1] 

(14) \\u-u\\ ,rX<^P-Hi,2,T for ueX{T), 
where 1 < r < 2 and u = mT 1 J n u(x)dx, m = m(Q). Then 

For our result on the long-time behavior, we introduce the discrete relative entropy 

> 0, k > 0. 



E[n k \n*\ = Y,™{K)n k K \og (^j 



Theorem 4 (Long-time behavior of finite volume solutions). Let the assumptions of The- 
orem^ hold and let (nx>, Sx>) be a solution to -(J7]) . Then for all k > 0, 



E[n k+1 \n*} + At(l - C*) | Vn^l^ + ^\\S k+1 " S*\\l,2,T < E[~ k 



n \n 



where C* = /i 2 C(fi) 2 ||no||L 1 (n)/(^£); C(Q) > only depending onfl, and ^ is the parameter 
in Q. In particular, if C* < 1, the discrete relative entropy is nonincreasing and 

-k/2 



\n k 



n*\\o,i,T < ^E[n°\n*}\\n \\ L i {n) (l + ^T=- At ) 



l^ll ttr < v '£^( 1 + I^ u )-™\ , eN , 
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where Cl > is the constant in the discrete logarithmic Sobolev inequality (see Proposition 
In particular, for each K G T, 

(n k K , S k K ) (n*, S*) as k ->• oo. 

In [22], the long-time behavior of solutions is shown under the condition /^H^olU 1 ^)/^ < 
C(Q) for a constant C(Q) > appearing in some Poincare-Sobolev inequality. We observe 
that our condition depends additionally on the regularity of the finite volume mesh. The 
convergence of \\n k — n*||o,i,r is approximately exponential for sufficiently small At > 0. 
The convergence result for \\S k — S*\\i t 2,T is weaker in view of the factor (At) -1 . In [22], 
the exponential convergence of S — S* is shown in the L 2 (Q) norm only and therefore, our 
result is not surprising. 

3. Existence of finite volume solutions 

In this section, we prove Theorem [T} The proof is based on the Brouwer fixed-point 
theorem. Let k G {1, . . . , M T — 1} and let (n fc , S k ) be a solution to ^ and (j7|), with k + 1 
replaced by k, satisfying (|9|)-([Io|). We introduce the set 

Z = {u G X(T) : u > in Q, ||u||z,i(n) < ||n ||i,i(n)}- 

The finite-dimensional space Z is convex and compact. In the following, we define the 
fixed-point operator by solving a linearized problem. First, we construct S G X(T) using 
the following scheme: 

(15) - J2 T - D SK,a + m(K)S K = 5 J2 r * Dn K,« + »m(K)n k K , KeT. 
Second, we compute n G X{T) using the scheme 

(16) ^^^k ~n k K )-J2 r*Dn Kja + ^ r a ((DS K>a ) + n K - (DS K , a )-n L ) = 0. 

a=K\L 



Step 1: Existence and uniqueness for (15) and (16). The linear system (15) can be 
written as AS = b, where A is the matrix with the elements 

A k ,k = r a + m(K), A K)L = -r a for K,L G T with a = K\L G £- mt ,K, 

a&£ K 

and b is the vector with the elements 

b K = S T ° Dn K,* + l^{K)n k K . 

aE£ K 

Since for all L G T, 

\A L , L \ ~ \ AK > L \ = m ( L ) > °> 

the matrix A is strictly diagonally dominant with respect to the columns and hence, A is 
invertible. This shows the unique solvability of (15). 
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Similarly, ( 16 ) can be written as Bn = c, where B is the matrix with the elements 
m(K) 



B 



K,K 



At 



Y t c {1 + {DS Kj(T ) + ), KeT, 

ae£ K 

B K ,L = -T a {l + (DSk,o)~) 



for K,LeT with a = K\L G & 



int,isr> 



and c is the vector with the elements ck = m.(K)n k K / At, K G T. The diagonal elements 
of B are positive, and the off-diagonal elements are nonpositive. Moreover, B is strictly 
diagonal dominant with respect to the columns since for all a = K\L G £\nt,Ki we have 
DSl )(T = —DSk,u which yields (DSl^) + = (DSk,v)~ and hence, 



\B 



L,L\ 



El* 



K,L\ 



m(L) 
At 



> 0. 



We infer that B is an M-matrix and invertible, which gives the existence and uniqueness 



of a solution n to (16) 



The M-matrix property of B implies that B~ x is positive. As a consequence, since n k and 
c are nonnegative componentwise, by the induction hypothesis, n = B~ x c is nonnegative 
componentwise. This means that n satisfies (j9|. Summing (16) over K G T, we compute 



Y m(K)n K = Y m(K)n k K = ||n |Ui(n)- 



KeT 



KeT 



Step 2: Continuity of the fixed-point operator. The solution to (15) and (16) defines 
the fixed-point operator F : Z — > Z, F(n 
For this, let 



n. We have to show that F is continuous. 



n 7 ) 7eN C Z be a sequence converging to n in X(T) as 7 — > 00. Setting 



77 ' 



F(n' y ) and n = F{n), we have to prove that n 1 — > n in X(T)- Using the scheme (15), 



we construct first S 1 (respectively, S) from n 7 (respectively, n). Then, using the scheme 
(16), we obtain n 7 (respectively, n). 
We claim that S 1 — S — > in X(T) as 7 — > 00. Indeed, since the map n \-t S, where 



S is constructed from (15), is linear on the finite dimensional space X(T), it is obviously 
continuous. Moreover, using the scheme (16) and performing the same computations as in 



the proof of Theorem 2.1 in [16], it follows that 



Y^{K)\nl-n K \ <2At(J2 



KeT 



\ n K\ 



KeT 



1/2 



t*\D(S^ - S) K , a \ 



KeT ae£ K 



1/2 



The right-hand side converges to zero as S* 7 — » S in X(T), which proves that n 7 — > n in 
•VCT). 

Step 3: Application of the fixed-point theorem. The assumptions of the Brouwer fixed- 
point theorem are satisfied, implying the existence of a fixed point of F, i.e. of a solution 
n k+1 to ^ satisfying (|9]). We have shown in Step 1 that (10) holds for n k+1 . Finally, we 
construct S k+1 using scheme (|7|. 
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4. A PRIORI ESTIMATES 

The proof of Theorem [2] is based on suitable a priori estimates which are shown in this 
section. We introduce a discrete version of the entropy functional used in [22] : 

E k = m ( K ) H ( n K), where H ( s ) = «(logs - 1) + 1. 

Proposition 5 (Entropy stability). There exists a constant C > only depending on Q, 
fJ', H^olU^n); an d £ (see Q ) such that for all k > 0, 



% -n K ) 



E k+1 ~E k + Y E ra\D(V^) K ^\ a + ™(K)\S k+1 \ 

(17) r,\DS^<CAt. 

Proof. By the convexity of H, we find that 

E k+i _ E k = J2 m (K)(#(r4 +1 ) - H(n k K )) < m ( K ) log(n^ +1 )(- fc+] 
KeT KeT 

Inserting the scheme (|6]), we can write 

E k+1 -E k <At r a (n k L +1 -n k K +1 )\ogn k K +1 

(T=K\L 

(18) - At J2 TriPS^n™ - (DS k + a 1 rn k L +1 ) logn^ 1 =: h + J 2 . 

<r=if|,L 

Now, we argue similarly as in the proof of Lemma 3.1 in [16J. We employ the symmetry 
of T a and a Taylor expansion of log around n k ^ 1 to infer that 



h = "At T ^ n K +1 ~ < +1 )(hgn k K +1 - \ogn k L 



cr=K\L 



ct=A"|L cr=X|L 

where fi^ +1 = t a n k ^ 1 + (1 — t^n^ 1 for some t CT G (0, 1). We perform a summation by parts 
in I 2 , using again the symmetry of r CT : 



h = -At T„((DS k £) + n k K +1 - (D^-n^Oogn^-logn 



L / 



c=K\L 
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Reordering the sum and using the expression for n k+1 in the Taylor expansion of log, it is 
shown in [TBI p. 468] that 



I 2 <-At r^DS^ilogn^ 1 -\ogn\ 



a=K\L 



The Taylor expansion shows that n^ +1 (logn^ 1 — logn^ +1 ) = n^ 1 — which gives 
I 2 <-At E r a DS k £(n k K +1 - n k L +1 ) = At E ^DS^Dn^. 

a=K\L a=K\L 



Summarizing the estimates for I\ and I2, (18) leads to 

(19) ^ - E k < -At £ r a + At £ r^^Dng. 

o-^IL <t= j FC|Z, 

The first term can be estimated for a = K\L as follows: 



fc+H Jr,!* 1 + J-ntf 1 



Drip* \ n L + V % 

(20) = V _^ |D(v^)^| > \D(V^) Kia \. 

In order to bound the second term, we multiply the scheme (|7j) by (At/8)S 1 ^~ l and sum 
over K e T: 

At 



By summation by parts, we find that 

At e ^^K = -f E ^Si 2 + TE m W< +1 ^ +1 

a=K\L a=K\L 

(21) -f E m wi^ +1 i 2 - 



Inserting (21) into (19) and employing (20), it follows that 

At 



E h+i_ E k + At r a \D(V^) K ^\ 2 + — E m Wis 



T&+1|2 

£ — J Jf 



(22) +f E < E 



a=K\L 
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It remains to estimate the right-hand side. We follow the proof of Proposition 3.1 in 
The Holder inequality yields 

(23) /i m W4 +1 4 +1 < /i||^ +1 ||o,6/5,rll^ +1 || ,6,r. 

KeT 

The discrete L 6 norm of S k+1 can be bounded by the discrete H 1 norm using the discrete 
Sobolev inequality (13) 

ll ak+l n / ofe+1 1| 

|l,2,Tj 



\S k+1 \\ 0fi ,T < C\\S k+1 \ 



where C > depends only on Q and £. For the discrete L 6//5 norm of n k+l , we employ the 
discrete Gagliardo-Nirenberg inequality [U Theorem 3] with 6 = 1/6: 



|^ +1 ||o, 6 /5,r = ||V^+ I II?^,^< C'||v / ^+ I ||?!^ ) ||\/^+ I 



I 20 

0,12/5,T — ~ II v 1 1 0,2,7" II v 1 1 1,2,7"' 



where C > depends only on f2 and £. Mass conservation (10) implies that 

||^ +1 ||o, 6 /5,r < C||n ||iT ( % (ikU^ + \^ r \ 2T T ■ 
With these estimates, (23) becomes 

29 



KeT 

< 2Cfi\\S k+ %^r\\n \\ 1 L - 1 d {n) (iKIU^) + | Vn^l^)' 
Then, by Young's inequality with p x = 2, p 2 = 2/(1 — 29), p 3 = 1/9, 



KeT 

+ l(\\n \\ L Hn) + \^ kTI \l 2 . 
< Jll^ 1 ||? Ar +||>/^|; ar + C'( < y ) Ai,||Tio|Ui( n )). 



Together with (22), this finishes the proof. □ 



Summing (17) over k = 0, . . . , My — 1 and using the mass conservation (10), we conclude 
immediately the following //-uniform bounds for the family of solutions (n v , S v ) v> o to (§- 
^ with discretizations T> v : 

(24) (n„), (n^logn,,) are bounded in L°°(0,T; L\tt)), 

(25) (V%) is bounded in L 2 (tt T ), 

(26) is bounded in L 2 (0,T; H\n)). 

Proposition 6. The family (n^) J)>0 is bounded in L 2 (0, T; if 1 (fi)). 
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Proof. First, we claim that (n v ) is bounded in L (0,T; W ' (f2)). To simplify the notation, 
we write n k+1 := n v (-,t k+1 ) G X(T n ). Applying the Cauchy-Schwarz inequality, we obtain 



K +1 ki,r„ = E m ^ n 



k+l ^k+l 



11 



K 



a=K\L 



< E E v^l^v^Wl ■ \/m{(T)d a 

1A t, ( E ( E ™(-k)4 +1 



n 



fc+i 
a: 



1/2 



Observe that in two space dimensions, 

m(cr)d(x^r, a) = 2m(K), 

rr££ K 

since the straight line between xk and xl is orthogonal to the edge a = K\L. Using this 
property, the mesh regularity assumption (Hj), and the mass conservation (10), it follows 
that 



1^ 1 1,1,7^ S 



(27) 



1/2 



1/2 



1/2 



n k+l 

v 1 1,2,7; 



E m (*>* 1 



^ +1 li,2,rJl n °lli 1 W 



In view of the entropy stability estimate (17), we infer that (n v ) is bounded in L 2 (0, T; 
W 1,1 ^)). Because of the discrete Sobolev inequality [H Theorem 4], 

n„fc+i 11 ^ ^ii/^^+i 11 

\\n v \\o,2,r v S c \\n v \\i,i,T n , 

the family (n v ) is bounded in L 2 (Q T ). 



In order to estimate the approximate gradient of n v , we employ (19). The last term in 
(19) is treated as follows. We multiply ([7]) by Atn^ 1 " 1 , sum over K e T, and sum by parts: 

At £ r^^ 1 ^ = -5At £ r.pn^i 2 



o-=i<:|L 



0"£^int, 
o-=K"|L 



+ //At £ m(K)|4 +1 | 2 - At £ m(K)S* +1 < +1 



<-5At T„\Dn k +l\ 2 + C{ 



n 



fc+l||2 

n IIo,2,t; 



, || afe+l||2 \ 
- n ^\\ a ri \\0,2,T n )- 



a=K\L 
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Inserting this estimate into (fl9|, we infer that 



E 



k+1 



E k 



At 

a=K\L 



r a \D(Vn^) K / + 5At £ T„\Dn k £\ 



CT=K"|L 



<^(ll< +1 Ho,2,r, 



I II Qfc+l||2 \ 

n ' \\ J V \\0,2,T V )- 



Summing this inequality over k = 0, . . . , Mt — 1 and observing that the right-hand side 
is uniformly bounded, we conclude that (Vn,) is bounded in L 2 (Qt), which finishes the 
proof. □ 

5. Convergence of the finite volume scheme 

We prove Theorem]!] Consider the family (n v , 5 , r? ) r? >o of approximate solutions to ([Sj)- (J7j) . 
In order to apply compactness results, we need to control the difference n v (-, t + r) — n^-, t). 
To this end, let E L°°(0,T; H 2+£ (tt)), where e > 0. We denote by § K the average of 6 
in the control volume K. Using scheme ^ and the notation n k+1 := n v (-,t k+1 '' 



At 



n 



,k+l 
K 



DS k + l a \)\D<p K A 



KeT v ae£ K 
< At(\n k+ \ 2 ,rM^, % 



fc+l|| lofc+ll 

1 0,2,-77, 1 11,2,^1^11,00,7;, 



,101 l,oo,7; ; ) 



■„ \i,2,r v + ||^1o,2,tJ^ +1 |i,2,tJ||0II# 2 +^), 



< CAt(|nJ +] - 11 " A - +1 i 
where C > only depends on Q. Summing over k — 0, . . . , Mt — 1 and employing Holder's 



inequality, the uniform bound on from Proposition [6] and on S v from (26) imply the 
existence of a constant C > 0, independent of r\, such that 

M T -1 

( 28 ) E E m(ir)(^ +1 -n^<L7At||0|| L 

°°(0,T;H 2 + £ (n))- 

fe=0 A'e7;, 

Now, similarly as in the proof of Lemma 10.6 in [TS], for all < r < At, 

T-T f 

(n v (x, t + r) — n v (x, t))(f)(x, t)dx dt 



PL 

Jo Jn 



Mr— 1 r-T-T 

< E / Xk(t,t + r)dt ™(KKn k K +1 -n k K )<b K , 

k=0 ^° KeT„ 



Xk(t,t 



where 

1 if kAt e (t,t + r], 
if k At £ it, t + r}. 

Inserting (28) into the above inequality and observing that 

-T-T 

Xk(t,t + r)dt < t < At, 
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we infer that 

/ / (nrj(x,t + r) - n v (x,t))(j)(x,t)dxdt < C||0||L°°(o,T;tf 2+e (n))- 
Jo Jn 

This gives a uniform estimate for the time translations of n v in L 1 (0, T; (H 2+e (Q))') . Since 
the embedding H l (Q) <^-> L p (fi) is compact for all 1 < p < oo in two space dimensions, we 
conclude from the discrete Aubin lemma [13] that there exists a subsequence of (n v ), not 
relabeled, such that, as r\ — > 0, 

n v — >■ n strongly in L 2 (0, T; L p (f2)), p < oo. 

Furthermore, since (V\) is bounded in L 2 (f2 T ), there exists y G L 2 (fi T ) such that 

Vfi, — 51 ?/ weakly in L 2 (fi<r). 

It is shown in the proof of Lemma 4.4 in [TO] that y = Vn in the sense of distributions. 
The bound of (S^) in L 2 (0,T; implies the existence of a subsequence, which is not 

relabeled, such that 

S v — 5, VS 1 ,, — ^ z weakly in L 2 (f2T)- 

Again, it follows that z = VS* in the sense of distributions. 

The limit rj — > in the scheme ([5])- ([7]) is performed exactly as in the proofs of Proposi- 
tions 4.2 and 4.3 in [16], using the above convergence results and the fact that (n^ViS 1 ,) 
converges weakly to nVS in L 1 (f2^). Compared to [TB], we have to pass to the limit also 
in the additional cross-diffusion term which does not give any difficulty since this term is 



linear in n v . This shows that (n, S) solves the weak formulation ( 11 )-( 12 ), finishing the 
proof. 

6. Proof of the discrete logarithmic Sobolev inequality 
The proof follows [121 Lemma 2.1]. Set 

v = ^ U ; U) e X(T), 

\\U - u\\ ,2,T 

where m = m(Q) and u = mT 1 J a udx. Then J Q vdx = and m -1 J Q v 2 dx = 1. Using 
Jensen's inequality for the (probability) measure m~ l v 2 dx, we find that for q > 2, 

— / v 2 \og(v 2 )dx = — - — / \og( v q ~ 2 )(rrr l v 2 dx) < — - — log) / v 9 ~ 2 (m~ 1 v 2 dx) 
mJn Q - 2 J n q - 2 \J a 

Q i / — 111 1 1 2 \ ^ Q ( — 111 1 1 2 i \ 

= ~q~^> l0g ^ M^t) < — ^ \Mo,q,T - 

because of logx < x — 1 for x > 0. With the definition of v, this inequality becomes 
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By the discrete Sobolev inequality ( 13 ), we infer that for 2 < q < 2dj (d—2) (and 2 < q < oo 
if d < 2) 



f. _. 2 (u — u) , 

/ (U — U) log -r. -J7777 dX < 

Jn m \\ u - u \\o.2,r 



u-u) 2 q C s (qf, 

u 



u.2. r 9 - s 

9 /^(g) 2 



1.2.T 



1 1 ||w - u||o, a ,r- 



9-2 V t 

Inequality (4.2.19) in [TF] (adapted to domains with general measure) shows that 



f (U 1i\ 2 

w^log — — j— j-j — 7777 dx < / (u — w) 2 log — — j-j-; rTj2 dx + 2\\u — m||o,2,t- 

n m \\ u \\o,2,T m W u ~ M llo,2,r 



Hence, with the discrete Poincare inequality (14), 



2i a <r 9 g g(g) 2 | 12 1 ^C5(9) 2 , . A - ]U] _ ||2 



9 („ ,^2 , ^5(9) 2 Cp(2) 2 , g-4^ ,. 2 . , , 2 



and Proposition [3] follows. 



7. Long-time behavior 



In this section, we prove Theorem HI Similarly as in the proof of Proposition^ (see (19) 



and (20)), we have 

E[n k+1 \n*] - E[n k \n*} = E ™(K)(H(n k + l ) - H(n k K )) 



KeT 



(29) < -At J2 r«\D(V^) K J + At J2 TaDSjgDn™ 



a=K\L a=K\L 



In view of the identity S* = fin*, we can formulate the scheme ([7]) for all K E T as 
= r„D(S k+1 - S*) K ,„ + 5 T«Dn k +l + m(K)(iu(n k K +1 - n*) - (S k { +1 - S*)). 

Multiplying this equation by (S^ 1 — S*)/5 and summing over K G T gives 
= -] E r a \D(S k+1 - S*) K „\ 2 - ^DS k £Dn k £ 

a=K\L a=K\L 
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Replacing the last term in (29) by the above equation and using the Cauchy-Schwarz and 
Young inequalities, it follows that 

E\n k+l \n*} -E[n k \n*} + At|vV+ T |' ^+ ^\\S k+1 -S*\\ 2 12T 

11,2,/ A it' 



fiAt 



"£™(K)(n k K +1 -n*)(S k K +1 -S*) 
KeT 



2^ \\ n n llo,2,r + 25 llo,2,r- 



The second term on the right-hand side can be absorbed by the corresponding expression 
on the left-hand side. For the first term, we employ the continuous embedding of BV(Q) 
into L 2 (Q) in dimension 2 and the definition of the seminorm | • 1 1,1,7- [1, Theorem 2]: 

||^ +1 -^||o,2,r<C(fi)|n fc+1 | lilir . 



Then, using inequality (27), we can estimate: 

(30) ||n fc+1 -n*||^ )T <^(Q)lno|Ui(n)|v^|' 2ir . 
Setting C* = ii 2 C(Q) 2 \\n Q \\ L i^/(5^), this yields 

E[n k+1 \n*] - E[n k \n*] + At(l - C*) \ 2 T + ^\\S k+l - S% 2jT ^ °- 

To proceed, we assume that C* < 1. With the discrete logarithmic Sobolev inequality 
(Proposition [3]) , 

£[n fc+ V] < Cilvn^l^T-, 

we infer that 

/ 1 _ n* \ At 

(31) E[n k+1 \n*] \l + ~^^) ~ E[n k \rf] + — \\S k+1 - S% 2>T < 0, 

and hence, 

(32) E[n k \n*} < ^1 + Atj E[n°\n*]. 

Then, by the Csiszar-Kullback inequality [8, Prop. 3.1] (this result is valid in bounded 
domains too), 

\\n k - n*\\ 2 ,i, T < M\n \\ LHn) E[n k \n*] < 4||n Q |U 1(n ) (l + ^r=~ At ) E[n°\n*). 
Going back to (31), we find that 

\\S k+1 ~ S%2,r < f^V] ^Xt{ 1 + E[n°\n% 
which concludes the proof. 
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8. Numerical experiments 

In this section, we investigate the numerical convergence rates and give some exam- 
ples illuminating the long-time behavior of the finite volume solutions to nonhomogeneous 
steady states. 

8.1. Numerical convergence rates. We compute first the spatial convergence rate of 
the numerical scheme. We consider the system Q on the square = (— |, I) 2 . The time 
step is chosen to be At = 10~ 8 , and the final time is t = 10~ 4 . The initial data is the 
Gaussian 

M ( (x - x ) 2 + (y - y ) 2 ' 
n {x,y) = ^^exp 



2n9 



29 



where 9 = 1CT 2 , M = ||rio||ii(n) — 67r, and xo = yo = 0.1. The model parameters are 
5 = 1CT 3 and \i = 1. We compute the numerical solution on a sequence of square meshes. 
The coarsest mesh is composed of 4 x 4 squares. The sequence of meshes is obtained by 
dividing successively the size of the squares by 4. Then, the finest grid is made of 256 x 256 
squares. The LP error at time t is given by 



eAa 



\nAx{-,t) - n cx (-,t)\\ LPin) , 



where n& x represents the approximation of the cell density computed from a mesh of size 
Ax and n ex is the "exact" solution computed from a mesh with 256 x 256 squares (and 
with At = 10~ 8 ). The numerical scheme is said to be of order m if for all sufficiently small 
Ax > 0, it holds that e& x < C(Ax) m for some constant C > 0. Figure [l] shows that the 
convergence rates in the L 1 , L 2 , and L°° norms are around one. As expected, the scheme 
is of first order. 



10' 



10 



10" 



- • - L 1 Error 

- • - L 2 Error 
-■- L inf Error 
X 


_ . n 

_ . m- " 




,m~~ _ . « 

^ • " " ^ ^ - ~ 



10" 



10" 



10" 



10" 



Figure 1. Spatial convergence orders in the L 1 , L 2 , and L°° norm. 
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8.2. Decay rates. According to Theorem [4| the solution to the Keller-Segel system con- 
verges to the homogeneous steady state if /x or 1/5 are sufficiently small. We will verify 
this property experimentally. To this end, let Q = 



■\, \) 2 and 



n (x,y) 



M 
2^9 



cxp 



x 2 + y 2 
26~ 



where 9 = and M = 5tt. We compute the approximate solution on a 32 x 32 Cartesian 
grid, and we choose At = 2 • 10~ 4 . In Figure [2j we depict the temporal evolution of the 
relative entropy E*(t k ) = E[n k \n*} in semi-logarithmic scale. In all cases shown, the 
convergence seems to be of exponential rate. The rate becomes larger for larger values of 



5 or smaller values of /x which is in agreement with estimate (32). In fact, the constant C* 
is proportional to /x 2 /5 (see Theorem E| and the rate improves if fi 2 /S is smaller. 



Relative entropy 



10" 



10 





o 


5 = 




• 
□ 
V 


5 = 0.01 ' 
5 = 0.05 
5 = 0.1 


▼ n ** t*. 






v 

v 






^ D □ - D 


*• n° 





0.2 



0.4 



0.6 



0.8 



(a) (j, = 1. 



W 10 



Relative entropy 




(b) <5 = 1(T 3 



Figure 2. Relative entropy i?[n fc |n*] versus time t k in semi-logarithmic scale 
for various values of 5 and //. 



As a numerical check, we computed the evolution of the relative entropies for different 
grid sizes iV and different time step sizes At. Figure [3] shows that the decay rate does not 
depend on the time step or the mesh considered. 



8.3. Nonsymmetric initial data on a square. In this subsection, we explore the be- 
havior of the solutions to ([T]) for different values of 5. We choose f2 = (— |, |) 2 with a 
64 x 64 Cartesian grid, /x = 1, and At = 2 ■ 10~ 5 . We consider two nonsymmetric initial 
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Relative entropy 



10" 



10" 




(a) At = 2.10" 



Relative entropy 



10" 



10" 




(b) N = 16. 



Figure 3. Relative entropy .E[n fc |n*] versus time t k in semi-logarithmic scale 
for various mesh and time step sizes. 



functions with mass 67r: 
(33) 

6tt 

(34) 



Air 
2^9 



exp 



(x - xp) 2 + (y- ypf 
29 

(x - xp) 2 + (y- y f 
29 



2tt 
2^9 



cxp 



(x - xi) 2 + (y - y x f 
29 



where 9 = 10 , xq = yo = 0.1, and X\ — y% — —0.2 (see Figure |4j). 

We consider first the case 5 = 0, which corresponds to the classical parabolic-elliptic 
Keller-Segel system. In this case, our finite volume scheme coincides with that of [16]. We 
recall that solutions to the classical parabolic-elliptic model blow up in finite time if the 
initial mass satisfies M > Air [2Z] (in the non-radial case). The numerical results at a time 
just before the numerical blow-up are presented in Figure |5j We observe the blow-up of 
the cell density in finite time, and the blow-up occurs at the boundary, as expected. More 
precisely, it occurs at that corner which is closest to the global maximum of the initial 
datum. 

Next, we choose 5 = 10~ 3 and 5 = 10~ 2 . According to Theorem [lj the numerical 
solution exists for all time. This behavior is confirmed in Figure [6j where we show the cell 
density at time t — 5. At this time, the solution is very close to the steady state which 
is nonhomogeneous. We observe a smoothing effect of the cross-diffusion parameter 5; the 
cell density maximum decreases with increasing values of 5. 



8.4. Symmetric initial data on a square. We consider, as in the previous subsection, 
the domain fl = (— |, |) 2 with a 64 x 64 Cartesian grid, \i = 1, and At = 10~ 5 . Here, we 
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(a) Initial datum no i- (b) Initial datum no 2- 

Figure 4. Initial cell densities. 




(a) Initial datum no,i, t = 1. (b) Initial datum no, 2, t = 0.6. 

Figure 5. Cell density computed from nonsymmetric initial data with M 
6ir and 5 = 0. 



consider the radially symmetric initial datum 

(35) n , 3 (x,,) = — expf _ 
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(a) Initial datum n ,i, 5 = 10 




(b) Initial datum n 2 , 6 = 10 




(c) Initial datum no,i, <5 = 10 2 . (d) Initial datum n ,2, S = 10 2 . 

Figure 6. Cell density computed at i = 5 from nonsymmetric initial data 
with M = 67r for different values of S. 

with M = 207r and = 1CT 2 . Since M > 8tt and the initial datum is radially symmetric, 
we expect that the solution to the classical Keller-Segel model (5 = 0) blows up in finite 
time (261 130] . Figure [7] shows that this is indeed the case, and blow-up occurs in the center 
of the domain. 

In contrast to the classical Keller-Segel system, when taking 5 = 10~ 3 , the cell density 
peak moves to a corner of the domain and converges to a nonhomogeneous steady state 
(see Figure [8]). The time evolution of the L°° norm of the cell density shows an interesting 



A FINITE VOLUME SCHEME FOR A KELLER-SEGEL MODEL 



23 



[l.6 




Figure 7. Cell density at time t = 0.05 computed from to the radially 
symmetric initial datum n Q 3 with M = 20n and 5 = 0. 

behavior (see Figure [9]). We observe two distinct levels. The first one is reached almost 
instantaneously. The L°° norm stays almost constant and the cell density seems to stabilize 
at an intermediate symmetric state (Figure |8^i). After some time, the L°° norm increases 
sharply and the cell density peak moves to the boundary (Figure |8}d). Then the solution 
stabilizes again (Figure |8p). We note that we obtain the same steady state when using a 
Gaussian centered at (— 10~ 3 , — 10~ 3 ). 




(a) t = 0.6. (b) t = 0.73. (c) t = 5. 



Figure 8. Cell density computed from the radially symmetric initial datum 
n , 3 with M = 20tt and 5 = 10" 3 . 
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x 10 



Hn(t)ll 







Figure 9. Time evolution of ||n fe ||L°°(n) computed from the radially sym- 
metric initial datum n ,3 with M = 207r and 5 = 10~ 3 . 



Nonsymmetric initial data on a rectangle. We consider the domain Q = (—1, 1) 
^, ^) and compute the approximate solutions on a 128 x 64 Cartesian grid with At = 5- 

1, and we choose the initial data n 0j i and n ,2, defined 
0, the solution blows up in finite time and the blow 

1CT 3 , the approximate 



8.5. 

x(- 

10 5 . The secretion rate is again /i = 
in (pl)-([34|) with mass M = 6vr. If 5 



up occurs in a corner as in the square domain (see Figure 10). If 5 



solutions converge to a non-homogeneous steady state (Figure 11). Interestingly, before 
moving to the corner, the solution evolving from the nonsymmetric initial datum uq^ shows 
some intermediate behavior; see Figure [TTp . 

8.6. Symmetric initial data on a rectangle. The domain is still the rectangle Q = 
(— 1, 1) x (—1, |), we take a 128 x 64 Cartesian grid, /i = 1, and At = 10~ 5 . We choose 
the initial datum rio,3, defined in (35), with M = 20tt. Clearly, the approximate solution 
to the classical Keller-Segel model 5 = blows up in finite time in the center (0,0) of the 
rectangle. When 5 = 10~ 3 , the cell density peak first moves to the closest boundary point 



before moving to a corner of the domain, as in the square domain (Figure 12). However 



in contrast to the case of a square domain, there exist two intermediate states, one up to 
time t « 0.9 and another in the interval (0.9,2.3), and one final state for long times (see 
Figure 13). We note that the same qualitative behavior is obtained using 5 = 10 -2 . 
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Figure 10. Cell density computed from nonsymmetric initial data with 
M = 6vr and 6 = 0. 
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